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1  Introduction 

This  is  the  final  technical  report  on  a  project  which  proposed  to  study  the  applicability 
of  several  time- frequency  representations  including  the  ambiguity  function  [7,  8,  9,  25, 
35],  Wigner-Ville  distribution  [33]  the  Gabor  transform  [19]  and  the  wavelet  transform 
[17,  18]  to  problems  in  nonlinear  filtering,  multichannel  signal  analysis,  and  detection 
and  feature  extraction  especial;  of  transient  signals  corrupted  by  non-Gaussian  noise. 
The  main  contributions  of  this  research  are 

•  The  development  of  adaptive  methods  for  selecting  windows  or  analyzing  signals 
targeted  to  a  given  application  which  lead  to  compact  representations  and  fast, 
numerically  stable  algorithms  for  analysis,  synthesis  and  processing. 

•  A  study  comparing  the  applicability  of  several  time-frequency  representations 
to  problems  in  signal  detection  and  feature  extraction  which  includes  the  effects 
of  such  representations  on  processing. 

•  The  design  and  comparison  of  code  implementing  time-frequency  signal  repre¬ 
sentation  and  synthesis  including  the  discovery  of  new  algorithms  with  signifi¬ 
cantly  different  behavior  from  those  in  standard  use, 

•  The  application  and  construction  of  optical  devices  for  implementing  and  pro¬ 
cessing  time-frequency  signal  expansions. 

Many  important  application  areas  including  radar  and  sonar,  seismology  and 
medical  imaging  require  tools  which  encode  and  process  joint  localized  time  and 
frequency  information.  Classical  Fourier  methods  average  information  taken  over  en¬ 
tire  signal  duration  and  as  such  do  not  directly  display  this  localized  information. 
Time-frequency  representations  offer  the  potential  of  analyzing  and  processing  signal 
information  in  a  space  where  this  localized  time-frequency  information  is  more  readily 
available  and  can  be  exploited  in  nonlinear  filtering,  pattern  recognition,  classification 
and  detection. 

Due  to  the  untimely  death  of  professor  George  Eichmann  of  CCNY,  professor 
Yao  Li  and  professor  Michael  Conner  of  CCNY  directed  most  of  the  research  into 
applications  with  direct  contract  support.  Professor  Izdor  Gertner  and  professor 
Joshua  Zeevi  of  Rutgers  University  played  collaborative  roles  in  formulating  specific 
applications  some  of  which  they  independently  tested.  The  mathematical  theory  as 
well  as  overall  direction  was  the  responsibility  of  professor  Richard  Tolimieri  of  CCNY 
who  received  additional  support  for  this  effort  from  DARPA  contract  ^F49620  89 
C0020.  Some  aspects  of  the  mathematical  research  were  carried  out  jointly  with 
Richard  Orr  of  Atlantic  Aerospace  Corporation. 
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2  Technical  Problem 

Classical  Fourier  methods  are  based  on  the  stationarity  assumption  on  the  signals  and 
Gaussian  assumption  on  noise.  Processing  techniques  are,  therefore  limited  in  most 
cases  to  the  application  of  time-invariant  systems  to  signals.  The  space  of  the  signal 
(the  time  domain)  and  the  space  of  its  Fourier  transform  (the  frequency  domain) 
contain  equivalent  information.  Filtering  can  be  implemented  in  either  domain  by 
convolving  the  signal  with  a  mask  in  the  time  domain  or  by  multiplying  the  Fourier 
transform  of  the  signal  by  a  filtering  function.  In  either  case,  the  effect  is  to  suppress, 
pass  through  or  enhance  certain  signal  frequencies  uniformly  across  the  signal.  Since 
the  Fourier  transform  averages  information  taken  over  all  time,  all  points  are  equally 
affected. 

The  stationarity  condition  is  not  valid  in  many  applications  and  does  not  faithfully 
represent  the  underlying  natural  process.  Traditional  techniques  based  on  Fourier 
analysis  do  not  exhibit  the  necessary  resolution  for  dealing  with  many  important 
signal  processing  tasks  including  those  required  in  sonar  processing  where  propaga¬ 
tion  and  multi-path  conditions  must  be  taken  into  account.  Time-frequency  repre¬ 
sentational  methods  provide  signal  representations  containing  joint  localized  time- 
frequency  information  offering  high  resolution  and  good  discrimination  capabilities 
when  non-stationarities  and  nonlinear  effects  distort  the  signal.  Processing  in  the 
space  of  this  representation  can  take  place  which  modifies,  selects,  and  distinguishes 
certain  frequencies  as  before  but  now  the  effects  of  this  processing  vary  over  time 
providing  powerful  tools  for  dealing  with  nonstationary  phenomena. 

A  major  distinction  between  Fourier  methods  and  time- frequency  methods  is  that 
in  the  latter  we  have  the  freedom  to  choose  a  localizing  window  or  analyzing  signal 
relative  to  which  signal  information  is  referenced.  The  choice  of  the  analyzing  signal 
can  affect 

•  The  existence,  uniqueness,  and  compactness  of  the  representation. 

#  The  speed  and  stability  of  computations. 

The  feasibility  of  time-frequency  methods  depends  to  a  large  extent  on  the  existence 
of  fast,  numerically  stable  algorithms  for  synthesis  and  analysis  of  signals  and  for 
processing  in  the  time-frequency  domain. 

As  in  the  classical  Fourier  analysis  of  signals,  performance  is  greatly  affected  by 
strategies  which  exploit  known  signal  information.  For  example  signal  information 
as  to  time  duration  and  bandwidth  are  essential  for  selecting  a  sampling  rate  which 
adequately  represents  the  signal.  The  matched  filter  is  the  optimum  signal  processing 
strategy  for  signal  detection  and  feature  extraction  in  noisy  environments  when  the 
signal  is  exactly  known.  Preprocessing  strategies  for  time-frequency  representations 
must  give  information  as  to  the  form  and  shaping  of  the  window  to  be  used  in  the 
analysis.  We  have  based  such  strategies  on  the  Zak  transform  and  the  ambiguity 
function. 
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3  Technical  Results 

We  have  attacked  the  problem  of  designing  efficient  time-frequency  computational 
tools  by 

•  Developing  selection  procedures  which  shape  an  analyzing  signal  from  a  priori 
and  precomputed  front-end  computations  on  input  data  based  on  Zak  transform 
and  ambiguity  function. 

•  Implementing  and  comparing  code  for  computing  Gabor  coefficients  based  on 
methods  found  in  [3,  11].  This  code  uses  fast  FFT  algorithms  developed  under 
DARPA  contract  #F49620  89  C0020.  We  have  determined  that  the  algorithm 
based  on  the  deconvolution  formula  in  [4]  produces  the  fastest  code  and  have 
applied  this  code  using  the  one-sided  exponential  window  to  transient  signal 
detection. 

•  Developed  a  new  algorithm  for  computing  classical  Gabor  coefficients  based 
on  the  concept  of  a  generalized  biorthogonal.  This  algorithm  ’’delays”  the  ef¬ 
fects  of  zero  theorem  and  provides  for  numerically  stable  computation  of  Gabor 
coefficients  locally  around  known  Gabor  coefficients. 

•  Developed  the  proper  form  of  finite  discrete  Gabor  transform  by  periodizing 
and  sampling  and  presented  the  results  in  [10,  32]. 

We  have  applied  these  results  to 

•  Gabor  representational  schemes  for  submicron  filtering,  image  reconstruction 
and  image  transfer  for  application  to  submicron  lithography 

•  Designing  and  constructing  optical  devices  to  implement  time- frequency  repre¬ 
sentations  and  to  carry  out  processing  on  such  representations  including  tran¬ 
sient  signal  detection. 

3.1  Zak  Transform  Selection  Procedures 

The  study  of  selection  procedures  in  the  case  of  thi  Gabor  transform  has  been  earned 
out  in  some  detail.  Any  Gabor  scheme  which  seeks  to  apply  to  all  square  integiable 
signals  will  be  limited  by  the  Heisenberg  uncertainty  principle.  The  zero  theorem 
[6]  and  more  generally  Dalian’s  theorem  are  concrete  ramifications  of  the  unceruainty 
principle  which  directly  obstruct  general  representation  schemes,  either  by  restricting 
the  class  of  analyzing  signals  or  by  requiring  oversampling  i.  e.  lattices  whose  fun¬ 
damental  region  area  is  less  than  one.  The  effort  in  this  research  was  directed  to  the 
design  of  adaptive  procedures  which  used  front-end  computation  of  the  Zak  trans¬ 
form  of  incoming  data  as  tools  for  choosing  optimal  windows.  The  Zak  transform  [36] 
has  played  a  major  role  in  this  research.  It  his  played  a  crucial  role  in  much  of  the 
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theory  and  application  of  time-frequency  representations  including  the  design  of  fast 
algorithms  for  determining  Gabor  coefficients  [3,  4, 11,  12].  The  numerical  stability  of 
these  algorithms  is  limited  by  the  Zero-theorem  [5].  These  algorithms  usually  involve 
division  by  the  Zak  transform  of  the  analyzing  signal  which  if  sufficiently  smooth 
must  have  a  zero.  Results  in  [1,  30]  suggests  methods  for  classifying  and  reducing  the 
negative  consequences  of  the  Zero-theorem. 

The  Zero-theorem  is  intimately  tied  to  the  Heisenberg  uncertainty  principle.  One 
of  the  first  proofs  [5]  rests  on  the  non-abelian  group  structure  of  the  Heisenberg  group. 
In  [21],  we  review  the  function  theory  on  the  Heisenberg  group,  more  accurately  on 
certain  compact  nil-manifolds  coming  from  the  Heisenberg  group  called  Heisenberg 
manifolds  and  apply  this  theory  to  the  problem  of  designing  adaptive  robust  Gabor 
procedures.  Many  of  these  results  are  special  to  the  Gaussian  and  make  use  of  detailed 
knowledge  of  the  Taylor  series  expansion  of  the  Zak  transform  of  the  Gaussian  at  its 
unique  zero  in  the  unit  square.  Starting  with  the  Gaussian,  we  can  define,  in  a 
natural  manner,  a  family  of  Gaussian-like  analyzing  signals  parameterized  by  the 
number  and  the  position  of  Zak  transform  zeroes.  In  fact,  several  families  can  be 
constructed  offering  some  variation  in  the  local  behavior  of  the  Zak  transform  at 
vanishing  points  and  global  envelope  characteristics. 

These  results  provide  powerful  tools  for  matching  analyzing  signal  to  application 
signal.  In  [5,  30],  this  led  to  a  detailed  description  of  the  oversampling  needed  to 
remove  the  numerical  instability  of  expanding  arbitrary  square  integrable  signals  as 
integer  time-frequency  translates  of  the  Gaussian.  A  more  exhaustive  analysis  [30] 
established  sufficient  conditions  on  signals  for  the  existence  of  good  Gabor  expansions. 

Recent  work  has  extended  these  results  to  signals  other  than  the  Gaussian  such 
as  Hermite  functions  [1],  one-sided  exponentials  and  Legendre  functions.  As  with  the 
Gaussian  a  family  of  analyzing  signals  can  be  associated  to  each  signal  of  a  type. 
The  analyzing  signals  in  a  family  can  as  before  be  parameterized  by  the  number  and 
the  position  of  Zak  transform  zeroes  however  as  we  move  through  the  families  we 
have  a  variety  of  local  behavior  at  the  Zak  transform  zeroes  and  global  Zak  transform 
envelope  characteristics. 


3.2  Ambiguity  Function  Selection  Procedures 

In  the  early  1960s,  research  summarized  in  [16]  into  the  Radar  ambiguity  function 
yielded  a  rather  a  complete  analysis  and  characterization  of  the  auto-  and  cross¬ 
ambiguity  function.  In  [24],  this  research  was  coupled  with  the  Poisson  summation 
formula  to  derive  fundamental  formulas  which  lie  at  the  heart  of  the  theory  of  Gabor 
expansions.  These  formulas  relate  the  energy  of  the  discrete  cross-ambiguity  function 
of  two  signals  over  a  lattice  with  the  inner  product  of  the  discrete  auto-ambiguity 
function  of  the  signals  over  a  complimentary  lattice.  The  lattice  in  the  former  defines 
the  lattice  of  time-frequency  translates  of  the  analyzing  signal  in  the  corresponding 
Gabor  expansion.  Frame  conditions  on  analyzing  signals  and  lattices  imply  a  good 
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representation  theory  for  all  square  integrable  functions.  We  will  extend  the  role 
played  by  these  fundamental  formulas  and  use  them  to  provide  tools  for  shaping 
good  analyzing  signals  and  determining  good  time-frequency  translation  lattices  for  a 
given  application.  Although  at  present  these  results  are  more  speculative  as  compared 
to  Zak  transform  results,  we  expect  that  significant  contributions  to  both  theory  and 
application  of  Gabor  expansions  will  be  achieved. 

3.3  New  Algorithms 

In  [11,  12,  13,  14],  the  Zak  transform  was  used  to  construct  a  biorthogonal  of  an 
analyzing  signal.  This  led  to  a  computational  procedure  for  computing  Gabor  co¬ 
efficients  by  taking  inner  products  relative  to  ai  biorthogonal  but  the  Zero  theorem 
once  again  created  problems.  For  analyzing  signals  having  continuous  Zak  trans¬ 
forms  the  biorthogonal  need  not  be  square-  integrable.  In  the  case  of  the  Gaussian, 
we  have  introduced  the  concept  of  a  generalized  biorthogonal  [32]  which  depending 
on  certain  moment  conditions  can  be  taken  as  smooth  as  desired  but  at  a  cost.  The 
smoother  the  generalized  biorthogonal  the  more  complex  and  numerically  unstable 
the  algorithm  computing  Gabor  coefficients.  Direct  inner  product  computations  are 
no  longer  sufficient  and  certain  difference  equations  must  solved.  However  Gabor  co¬ 
efficient  computation  is  localized  in  the  sense  that  if  a  Gabor  coefficient  is  known  then 
in  some  region  in  the  lattice  around  this  Gabor  coefficient,  other  Gabor  coefficients 
can  be  efficiently  computed.  The  extent  to  which  this  statement  can  be  made  precise 
and  the  role  of  such  an  algorithm  in  applications  is  currently  under  study. 

3.4  Optical  Processing 

One  aspect  of  this  program  was  trying  to  identify  the  role  of  optics  in  the  implemen¬ 
tation  and  processing  of  wavelets  and  other  promising  decompositions.  Through  our 
two  years  of  theoretical  and  experimental  investigations,  it  has  become  quite  clear 
that  for  many  signal  processing  and  radar  applications,  optics  may  play  an  important 
role  to  speed  up  the  entire  process.  Our  findings  have  been  published  in  five  refer¬ 
eed  journal  articles  with  one  additional  article  in  the  process  of  publication.  A  brief 
summary  of  the  findings  is  give  below. 

Gabor  expansion  of  one  dimensional  short  signal  sequences  has  been  generated 
optically  for  the  first  time.  Through  the  use  of  the  biorthogonal  signal  of  the  selected 
Gabor  window,  the  expansion  of  an  arbitrary  signal  upon  the  window  was  performed 
optically  using  a  modified  optical  ambiguity  function  generator.  In  our  experimental 
effort,  the  single-sided  exponential  function  was  selected  as  the  Gabor  window  for  the 
detection  of  transients.  An  acousto-optic  processor  to  implement  the  approach  was 
proposed  and  conceptually  demonstrated. 

Wavelet  expansion  which  can  overcome  many  disadvantages  inherent  to  the  Ga¬ 
bor  expansion  was  also  studied  for  its  optical  implementation  for  the  first  time.  It 
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was  found  that  optics  is  quite  suitable  to  generate  and  display  both  the  direct  and 
the  inverse  wavelet  transforms  in  parallel.  Unlike  the  digital  computer-based  imple¬ 
mentation  which  prefers  the  course- to- fine  serial  generation  procedure,  using  optics, 
the  one-dimensional  to  two-dimensional  coordinate  expansion  could  be  performed  in 
parallel  and  on-the-fly.  An  optical  processor  to  process  one  dimensional  short  data 
sequence  was  first  built  and  tested.  The  processor  contains  one  DC  band  and  four 
other  higher  frequency  channels.  The  decomposition  of  a  spatial  chirp  signal  into 
these  bands  and  the  reconstruction  of  the  chirp  using  the  wavelet  decomposed  signals 
were  successfully  performed  for  the  first  time. 

A  further  step  was  taken  toward  identifying  the  suitability  of  using  optics  for  the 
multichannel  signal  analysis.  Both  the  Gabor  and  the  wavelet  transforms  were  studied 
in  terms  of  their  complexity  of  optical  implementations  in  general.  Both  one  and  two 
dimensional  signal  cases  were  assumed.  It  was  found  that  the  role  of  optics  in  the 
Gabor  case  was  to  process  the  ambiguity  function  necessary  for  the  implementation. 
While  for  the  signal  recovery,  although  it  is  not  impossible,  it  is  very  difficult  to  use 
the  optical  method.  Also,  due  to  the  use  of  coherent  processing,  it  is  the  square  of, 
not  the  Gabor  coefficients  themselves,  which  is  finally  generated.  Thus,  the  optical 
methods  are  more  suitable  for  the  detection  of  certain  signal  signatures  based  on 
the  selected  window  function  than  using  it  as  a  mathematical  expansion  tool.  On 
the  other  hand,  optics  finds  itself  promising  for  generating  both  direct  and  inverse 
wavelet  transforms.  This  is  the  case  for  processing  both  the  one  and  two  dimensional 
signals.  A  detail  comparison  of  the  space- band  width  product  related,  filter  dynamic 
range  limited,  as  well  as  the  misalignment  caused  performance  degradations  for  both 
applications  was  performed. 

To  extend  the  method  of  processing  the  Gabor  and  wavelet  transforms  of  one 
dimensional  short  signal  sequence  to  a  more  practical  situation  where  the  one  dimen¬ 
sional  long  data  sequence  are  encountered,  a  novel  optical  processor  concept  which 
scans  the  input  long  data  sequence  into  a  two  dimensional  format  was  proposed  and 
numerically  simulated.  Using  the  Chinese  remainder  theory  with  the  unity  difference 
between  the  two  selected  relative  moduh,  the  scanning  pattern  can  be  physically  im¬ 
plemented  through  a  modification  of  any  raster-scan  display  device.  The  mapped 
signal  can  then  be  optically  processed  in  a  standard  optical  image  processor  for  its 
wavelet  or  Gabor  processing.  Advantages  and  shortcomings  were  analyzed. 

3.5  Applications  to  Submicron  Lithography 

Analytic  and  algorithmic  methodologies  were  developed  for  submicron  filtering,  image 
reconstruction,  and  image  transfer  in  vector  and  scalar  forms  for  study  as  to  their 
applicability  to  micron  lithography.  This  study  contained  two  parts.  In  the  first  a 
projected  signal  had  to  be  addressed  in  a  scalar  as  well  as  a  vector  form  and  a  Gabor 
representation  of  the  final  line  profile  had  to  be  computed.  The  second  task  included 
a  representation  of  electron  beam  applied  to  submicron  masks,  such  that  optimal 
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intensity  evaluation  at  each  pixel  will  allow  analysis  of  improved  resolution. 

We  will  describe  in  detail  the  several  stages  of  the  study. 

In  first  part,  we  developed  and  formulated  the  algorithm  or  the  aerial  image 
methodology  in  a  scalar  form  for  partially  coherent  light  with  .5  micron  line  width. 
The  partial  coherence  was  addressed  through  the  ’’mutual  coherence  function”  that 
correlates  the  various  light  source  elements.  After  performing  a  stationary  phase 
calculation,  we  obtained  this  function  in  terms  of  Bessel  function  of  order  1.  Each 
entry  in  the  four-fold  integral  final  result  is  a  two-fold  integral.  Thus,  a  fast  algorithm 
was  needed  to  avoid  extremely  expensive  computations.  Variable  Gauss-Legendre 
quadrature  with  accuracy  control  coupled  with  the  interior  integrations  order  was 
chosen  and  implemented. 

In  the  next  stage,  the  aerial  image  for  vectorial  electric  field  was  investigated  and 
implemented.  The  rays  were  traced  through  the  entrance  as  well  as  the  exit  pupil. 
The  system  is  more  complicated  due  to  two  extra  integrations.  We  decided  to  lower 
the  order  of  the  innermost  loops  and  check  their  accuracy  via  an  averaged  monte-carlo 
simulation  of  the  integral.  The  polarization  partial  coherence  matrix  was  formulated 
and  the  necessary  algorithm  developed  and  tested. 

In  the  next  stage,  attention  was  given  to  the  exposure  systems  of  electron  beams. 
We  expressed  the  aerial  image  of  E-beam  as  a  linear  combination  of  double  Gaussian 
in  each  pixel  due  to  back-scattering.  The  aerial  image  is  thus  expressed  as  a  Gabor 
expansion  and  is  ready  made  for  data  compression  algorithm.  The  Haar  basis  and  its 
corresponding  wavelet  basis  failed  to  achieve  optimal  response.  A  projection  system 
has  been  implemented,  using  smooth  ”hat”  functions  basis  with  correct  polarity,  thus 
allowing  optimal  resolution  and  data  compression.  A  new  data  compression  algorithm 
was  developed. 

The  aerial  image  algorithms  were  generalized  to  include  models  with  time-dependent 
absorption.  The  mechanism  employed  is  a  time  dependent  filter.  In  other  words,  the 
filter  is  being  degraded  gradually  within  its  lifetime.  The  kinetics  of  the  filter  is  ob¬ 
tained  via  fifth  order  Adams-Bashforth  algorithm,  and  the  field  is  refiltered  at  each 
time  step.  The  necessary  algorithm  has  been  developed  and  implemented. 

A  vector  Maxwell  equation  solver  has  been  developed  using  spectral  elements 
methods.  An  image  is  constructed  and  is  deblumed  through  the  dynamic  filter  devel¬ 
oped  earlier.  The  moving  format  is  expressed  as  an  eikonal  system  with  a  curvature 
controlled  velocity  of  propagation,  and  this  formalism  lends  itself  to  self  improvement 
via  a  WKB  scheme  with  the  leading  term  obtained  from  the  eikonal  system.  The 
corresponding  PDEs  have  been  analyzed. 

The  previous  results  have  been  integrated  into  a  system  of  nonlinear  reaction- 
diffusion  equations  for  the  post-exposure  baking  and  a  second  system  of  dissolution. 
The  reaction  part  has  been  solved  analytically,  and  all  the  filtered  concentrations  are 
expressed  explicitly  in  terms  of  the  first  concentration,  which  in  turn  is  expressed 
implicitly.  A  very  fast  LV-decomp>osition  algorithm  has  been  developed  to  address 
finite-size  effects  while  diffusion  is  deconvolved  using  our  earlier  bases. 
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4  Summary  of  Results 

^  Code  in  Fortran  and  C  implementing  adaptive  selection  of  time-frequency  shifted 
Gaussian  and  Hermite  functions  as  analyzing  signals  relative  to  an  incoming  sig¬ 
nal. 

•  Code  testing  the  reliability  of  the  adaptive  selection  procedure  with  comparisons 
to  direct  standard  methods. 

•  Code  computing  generalized  biorthogonals  and  testing  the  numerical  stability 
and  robustness  of  Gabor  coefficient  computation. 

•  An  acousto-optic  processor  to  implement  Gabor  coefficient  computation  of  one¬ 
dimensional  short  signal  sequences.  The  single-sided  exponential  was  selected 
as  th'  Gabor  window  for  the  detection  of  transients. 

•  An  optical  processor  was  built  and  tested  to  compute  the  wavelet  expansion 
and  its  inverse. 

•  detail  comparison  of  Gabor  and  wavelet  transforms  was  undertaken  in  term  of 
their  complexity  of  optical  implementations. 

•  Methods  proposed  and  analyzed  for  optical  processing  of  Gabor  and  wavelet 
transforms  of  one-dimensional  data  sequences  using  the  Chinese  remainder  the¬ 
orem. 

•  Analytic  and  algorithmic  methodologies  were  developed  for  submicron  filter¬ 
ing,  image  reconstruction,  an  image  transfer  and  were  applied  to  submicron 
lithography. 

5  Announced  Results 

The  participants  in  this  project  have  published  their  results  in  many  papers  both  in 
refereed  journals  and  in  invited  articles  in  proceedings.  We  have  listed  most  significant 
below. 


5.1  Digital  Computations 

1.  L.Auslander,  M.An,  M. Conner  I.Gertner  and  R.Tolimieri,  Fine  Structure  of  the 
Classical  Gabor  Approximation,  IMA,  39  (2),  11-21,  Springer- Verlag,  1992. 


2.  L.Auslander,  C.Buffalano,  R.Orr  and  R.Tolimieri,  A  Comparison  of  the  Gabor 
and  Short-time  Fourier  Transforms  for  Signal  Detection  and  Feature  Extraction 
in  Noisy  Environments,  preprint,  1990. 


9 


3.  L.Auslander,  I.Gertner  and  R.Tolimieri,  The  Discrete  Zak  Transform  applica¬ 
tion  to  Time-frequency  Analysis  and  Synthesis  of  non-stationary  signals,  IEEE 
Trans.  Signal  Process.,  39  (4),  825-835,  April,  1991. 

4.  — ,  The  Finite  Zak  Transform  and  the  Finite  Fourier  Transform,  IMA,  39  (2), 
Springer- Verlag,  1992,  21-36. 

5.  L.Auslander,  and  R.Tolimieri,  On  Finite  Gabor  Expansions  of  Signals,  IMA 
Proceedings  on  Signal  Processing,  22,  Springer- Verlag,  Berlin/New  York,  1990. 

6.  Conner,M.,  Tolimieri,  R.,  and  Orr,  R.S.,  A  New  Algorithm  for  Computing  Gabor 
Coefficients  SPIE  San  Diego,  1991,  Special  edition  of  the  Journal  of  Optical 
Communications,  to  be  Published. 

7.  Orr,  R.S.  and  Tolimieri,  T.,  Poisson  Summation,  the  Ambiguity  Function  and 
the  Theory  of  Weyl-Heisenberg  Frames,  Submitted  to  IEEE  Trans.  IT,  July 
1990. 

8.  Tolimieri,  R.,  Problems  in  Gabor  Representation,  Proceedings  of  NATO,  ASI 
1992. 

5.2  Optical  Results 

1.  Li,  Y.,  Optics  for  Wavelet  Processing,  (Invited  Talk:  FA2)  Optical  Society  of 
America  Annual  Meeting,  Albuquerque,  NM,  Sept.  21-  25,  1992. 

2.  Zhang,  Y  and  Li,  Y.,  Optical  determination  of  Gabor  coefficients  of  transient 
signals.  Optics  Letters,  13  1991,  1031-1033. 

3.  — ,  An  optically  implementable  algorithm  for  convolution  correlation  of  long 
data  streams.  Optics  Communications,  85,  1991,  473-480. 

4.  Zhang,  Y.,  Li,  Y.,  Kanterakis,  E.  G.,  Katz,  A.,  Lu,  X.  J.  and  Tolimieri,  R. 
Optical  realization  of  wavelet  transform  for  a  one-  dimensional  signal.  Optics 
Letters,  14,  1992,  210-212. 

5.  Lu,  X.  J.,  Katz,  A.,  Kanterakis,  E.  G.,  Li,  Y.,  Zhang,  Y.,  and  Caviris,  N.  P. 
Image  analysis  via  optical  wavelet  transform.  Optics  Communications,  1992  in 
print. 

6.  Zhang,  Y.  and  Li,  Y.,  Coherent  optical  processing  of  Gabor  and  Wavelet  expan¬ 
sions  of  ID  and  2D  signals.  Optical  Engineering,  31,  1992,  1865-1885. 

7.  Zhang,  Y.  Li,  Y.  and  Caulfield,  H.  J.  Optical  implementation  of  long  data  stream 
convolution,  submitted  to  Applied  Optics  for  publication. 
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Appendix 

A  Zak  transform 

For  simplicity,  we  will  describe  several  basic  properties  of  the  Zak  transform  of  one¬ 
dimensional  signals.  Multidimensional  extensions  follow  by  the  same  arguments  [21, 
31,  34,  37,  38]. 

TheZaA:  trans/orm  of  a  signal  /  E  is  the  function  of  two  variables  Z(/)(x,y), 

X,  p  €  defined  by  the  formula 

k 

Z{f){x,p)  can  be  interpreted  as  the  discrete  Fourier  transform  at  x  of  the  sequence 
f{y  +  k),  k  E  Z  and  hence  contains  joint  time-frequency  information  about  the  signal 
/.  The  functional  equations 

^(/)  =  (a;  +  l,y)  =  Z{f){x,y) 

^(/)  =  (x,y  +  l)  =  e-2-^Z(/)(x,j/) 

imply  that  the  Zak  transform  is  completely  determined  by  its  values  on  the  unit 
square.  We  define  the  inner  product  of  two  Zak  transforms  by  the  formula 

<  Z(/i),Z(/2)  >=  /'  ['  Z{f^){x,y)Z*if,){x,y)dxdy. 

Jo  Jo 

Theorem  1  For  fi  and  f2  in 

</l,/2>  =  <  Z(/,),Z(/2)> 

=  EE  j[,‘  Mv+Miy  +  f 

j  k  ° 

=  S/  My  +  j)f2iy  +  k)dy 

j 

=  /  My)f2iy)(^y 

J—oo 

=  </l,/2>. 

The  following  theorem  characterizes  the  space  of  Zak  transforms. 


Proof: 

<  Z(/i).  Z(f2)  > 
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Theorem  2  If 

[  f  \  F{x,y)  f  dx4y  ^  or. 

Jo  Jo 

F{x  +  l,y)  =  F{x,y) 

Fix,y-\-l)  =  e-^^'^Fix,y) 
then  F  is  the  Zak  Fansform  of  a  unique  function  f  € 

Proof:  By  the  functional  equations,  we  have 

r 

where 

/r(«^4  1)  =  /r.fl(2/)- 

The  Zak  transform  and  the  Fourier  transform  are  related  by  the  following  formula. 

Theorem  3 

Proof:  Apply  the  proceeding  formula  to 

a{x,y)  =  e-^’‘'-‘Z{f)l-y,x). 

Define 

fmn{t)  =  fit-n)e^^'^\  m,n€Z. 

Theorem  4 


Proof: 


2-nikx 


=  02  f{y  + 

k 

The  application  of  the  Zak  transform  to  time-frequency  representation  is  the  con¬ 
sequence  of  two  fundamental  formulas  which  we  state  without  proof. 

1st  Fundamental  Formula 


Z(f){x,y)Z-U')(x,y)  =  E  E  <  /.S™.  > 

m€Z  n€Z 


Applications 
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•  Z{f){x,y)  vanishes  only  on  a  set  of  measure  0  implies  the  set 

{5mn  :  m,  n  e  Z} 


spans  Z/^(0J). 

•  I  Z{g){x,y)  1=  1,  I'.e.,  iff  the  set  {gmn  '■  m,  n  G  Z}  is  an  orthonormal  basis  of 

The  second  fundamental  formula  is  directly  tied  to  Gabor  expansions,  i.  e.,  ex¬ 
pansions  of  the  form 

f  ~  ^mnQmn'  (1) 

meZneZ 

We  will  assume  that  the  coefficients  in  (1)  satisfy 

51  2  I  imn  r<  OO-  (2) 

m€Z  n€Z 

In  this  case,  the  series  in  (1)  converges  to  /  in  Z/^(3i). 

Applying  the  Zak  transform  to  both  sides  of  (1)  we  have 

2nd  Fundamental  Formula 


Z(/)(x,!()  =  Z(a)(x,!,)  x:  x: 

m€Z  n^Z 

The  series  converges  in  L^(P)  to  a  doubly  periodic  norm-square  summable  function 
over  the  unit  square.  Multiplying  both  sides  by  Z*{g){x,y),  we  have  the  following. 

EE  <  /.fc"  >  =1  Z(g)(x,y)  I"  EE”-'"'''”*'"’"’'’-  (3) 

m  n  m  n 

The  right  hand  side  can  be  viewed  as  a  convolution. 

For  a  fixed  analyzing  signal  g,  a  signal  /  €  L^{dt)  need  not  have  an  expansion 
satisfying  (2).  The  second  fundamental  formula  is  the  main  tool  for  deriving  criteria 
on  /  and  g  guaranteeing  such  an  expansion  and  for  designing  algorithms  computing 
the  Gabor  coefficients  m,n  €  Z.  In  fact,  if  the  quotient 


is  norm-square  summable  over  the  unit  square  then  the  Fourier  coefficients 

bmn  =  I  I  (5) 

JlJl  Z{g){x,y) 
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satisfy  (1)  and  (2).  The  problem  is  that  if  Z{g)  is  continuous  then  it  must  vanish 
at  some  point,  a  result  proved  in  [5]  and  henceforth  called  the  Zero  theorem.  The 
quotient  (4)  is  not  necessarily  norm-square  summable. 

The  Zero  theorem  also  lies  at  the  heart  of  the  synthesis  problem  of  a  signal  from 
its  short-time  Fourier  transform.  The  inner  products,  <  /,5'mn  >,  m,  n  €  Z,  de¬ 
termine  the  product  Z{f)Z*{g)  a.  e.  but  division  by  Z{g)  t'  recover  Z{f)  is  not 
necessarily  numerically  stable.  For  example,  uniformly  sampluig  both  sides  of  the 
first  fundamental  formula  over  the  unit  square  can  produce  samples  of  the  product 
in  terms  of  the  2-dimensional  finite  Fourier  transform  of  periodized  inner  products 
but  as  we  increase  the  resolution,  zeros  of  Z{g)  will  be  approached  whenever  Z{g)  is 
continuous.  Nonuniform  sam’.'>ling  to  avoid  these  zeros  results  in  numerically  unstable 
computation  of  product  samples. 

For  fixed  g  G  with  Z{g)  continuous,  there  exists  /  G  which  do  not 

admit  Gabor  expansions  of  the  form  (1)  and  which  cannot  be  resynthesized  from 
short-time  Fourier  transform.  The  special  case  of  ^  was  considered  in  detail 

in  [5,  30],  where  it  was  found  that  every  /  G  with  Z{f)  continuous  has  a 

Gabor  expansion  if  we  allow  half  integer  shifts  in  the  time  or  frequency  variable,  i.e., 
n  ranges  over  ^Z,  The  main  results  in  [30]  give  a  precise  count  of  the  number  of 
half  integer  shifts  required  to  maintain  various  degrees  of  smoothness  in  the  Gabor 
expansion.  In  [18],  along  with  many  other  results,  an  account  of  similar  ideas  were 
presented  in  the  language  of  frames. 

In  the  following  subsections,  we  will  review  more  deeply  these  ideas  and  extend 
them  to  provide  tools  for  an  adaptive  Gabor  theory  and  for  the  design  of  efficient 
Gabor  expansion  implementation. 

A.l  Analyzing  Signal  Parameters 

Two  important  analyzing  signal  parameters  are 

•  the  zero  set  of  their  Zak  transform 

i*  the  deviation  of  the  absolute  value  of  their  Zak  transforms  from  unity. 

The  zero  set  affects  the  existence  and  compactness  of  Gabor  expansions  while  the 
deviation  is  a  measure  of  the  defect  from  orthonormality  of  the  wavelet  system. 

The  results  in  [5]  provide  a  parameterized  family  of  analyzing  signals  characterized 
by  having  Zak  transforms  with  unique  ‘analytic’  zeros  in  the  unit  square.  These 
results  can  be  extended  to  define  larger  families  of  Gaussian-based  analyzing  signals 
having  more  intricate  and  pre-assigned  zero  sets  and  indicate  general  methods  for 
constructing  non-Gaussian  based  analyzing  signals. 

The  Zak  transform  of  the  Gaussian  g{i)  =  can  be  written  as 


Zi9){x,y)  =  e 
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where  v  is  the  classical  theta  function  of  characteristic  (0,0) 

=  E  z  =  x  +  iy. 

rez 

It  is  well-known  that  v  is  an  entire  function  having  a  unique  zero  in  the  unit  square 
at  re  =  y  =  |. 

In  general,  if  guv{t)  =  g{t  —  u,  u  €  3?,  then 

Z{guv){x,  y)  =  -f  u,  y  -  u).  (6) 

Applying  (  6)  to  the  Gaussian  g{t)  =  we  see  ihat  the  collection 

:  0<u,u<l} 

of  time-frequency  translates  of  the  Gaussian  have  Zak  transforms  with  unique  ‘ana¬ 
lytic  zeroes’  in  the  unit  square  and  each  point  in  the  unit  square  is  the  Zak  transform 
zero  of  exactly  one  function  in  the  collection. 

Gaussian-based  signals  having  more  complicated  zero  sets  can  be  constructed 
using  the  following  formula 

Z{gmn){x,  y)  =  y),  m,  n  €  Z.  (7) 

If  P(x,  y)  is  the  trigonometric  polynomial 

=  EE 

m  n 

where  the  double  summation  is  finite,  then  (  7)  implies 

P(x,y)Z(y)(x,y)  =  Z(X]E  ^mnd'mn  )(^,y).  (8) 

m  n 

The  zero  set  of  the  Zak  transform  of 
A  =  EE  ^mndrnn 

m  n 

is  the  union  of  the  zero  set  of  Z{g){x,y)  and  the  zero  set  of  the  trigonometric  poly¬ 
nomial  P{x,  y). 

This  result  can  be  used  to  construct  signals  having  preassigned  Zak  transform 
zero  sets.  For  fixed  g  with  Z(y)  continuous,  the  zero  set  of  Z{g)  is  also  included. 
Applied  to  the  Gaussian  case,  g{t)  =  we  can  build  an  extensive  collection  of 

signals  with  a  wide  variety  of  Zak  transform  zero  sets. 

By  creating  a  families  of  signals  having  a  wide  range  of  Zak  transform  zero  sets,  we 
introduce  into  the  selection  procedure  a  criteria  based  on  matching  signal  to  analyzing 
signal  Zak  transform  zero  sets.  One  measure  of  the  effectiveness  of  this  approach  will 
be  the  compactness  of  the  representation  and  the  efficiency  of  the  computation. 
Example  1  The  trigonometric  polynomial 
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P^{x,  j/)  =  1  + 

has  exactly  two  zeros  in  the  unit  square,  at  (0, 0)  and  ( j).  The  product  P\{x,  y)Z{g){x,  ?/), 
g  =  ,  has  a  second  order  zero  at  (j,  j)  and  a  first  order  zero  at  (0,0).  Translating 

g  gives  rise  to  a  product  which  has  three  first  order  zeros. 

Example  2  The  trigonometric  polynomial 

P2{x,y)  =  2  +  ie^^'^  +  ie'^^'y 

has  exactly  one  zero  in  the  unit  square,  at  (j,  j).  The  product  P2{x,y)Z{g){x,y), 
g  =  has  first  order  zeros  at  (j,  j)  and  (j,  |). 

The  general  results  (  6),  (  7)  and  (  8)  can  be  applied  to  signals  other  than  the 
Gaussian.  In  preliminary  studies,  we  have  carried  out  these  constructions  on  such 
‘naturally’  defined  signals  as  one-sided  exponentials,  Hermite  functions  and  finite 
interval  restrictions  of  periodic  signals.  We  intend  to  increase  the  list  of  naturally 
defined  signals  to  other  special  functions  and  to  include  digital  signals  which  have 
played  important  roles  in  applications. 

In  [5],  partition  of  unity  arguments  matched  to  the  underlying  Heisenberg  group 
structure  were  used  to  prove  several  important  results.  These  ideas  can  also  be  used  to 
construct  signals  whose  Zak  transforms  have  pre-assigned  zero  sets  and  pre-assigned 
Taylor  expansions.  The  feasibility  of  constructing  signals  in  this  way  is  under  study. 

In  general,  the  collection 

{^mn  ••  m,  n  €  Z} 

is  not  orthogonal,  as  one  can  see  from  the  Gaussian  g{t)  =  e""'*.  The  condition  for 
orthonormality 

\  Z{g){x,y)\=l,  a.e.,  (9) 

is  a  simple  consequence  of  the  first  fundamental  formula  applied  to  f  =  g: 

I  Z(g)(,x,y)  P=  EE  < 

m  n 

By  the  zero  theorem,  if  Z{g)  is  continuous  then  (  9)  can  not  hold.  A  more  exact 
statement  about  the  continuity  of  Z{g)  expressed  in  terms  of  rate  of  decay  of  g  and 
orthonormality  is  given  by  Dalian’s  theorem.  There  are  several  ways  of  measuring 
defect  from  orthonormality.  For  example 

/  /  \\Z{g){x,y)\^ -l\'^  dxdy  =  '^'^\<  g,g,nn>\'^, 

where  J2m  J2n  denotes  summation  over  all  m,n  E  Z  except  m  =  n  =  Q,  measures  the 
energy  in  the  norm-square  sense  of  the  inner  product  <  g,gmn  >,  m,n  €  Z,  except 
m  =  n  =  0. 
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A. 2  Selection  Procedure 

A  signal  /  has  a  Gabor  expansion 

f  =  Yl<^mngmn  (10) 

m^Z n£Z 

satisfying 

53  53  I  a*nn  P<  OO 

meZ  neZ 

if  and  only  if 

The  Gabor  coefficients  in  (  10)  are  the  Fourier  coefficients  of  (  11).  In  general,  the 
smoother  the  quotient,  the  more  rapid  the  decay  of  the  Gabor  coefficients  and  the 
more  ‘compact’  the  Gabor  expansion  (  10)  in  the  sense  that  finite  partial  sums  should 
better  represent  /.  In  [1],  a  computer  experiment  verified  these  results  for  the  special 
case  of  the  Gaussian  as  an  analyzing  signal.  In  fact,  it  was  shown  that  if  the  quotient 
is  not  in  then  a  ‘good’  Gabor  expansion  does  not  exist. 

Selection  Procedure  For  a  given  signal  /,  compute  the  Zak  transform  of  / 
and  determine  its  zero  set.  From  a  library,  choose  g's  whose  Zak  transform  zero  set 
is  contained  in  the  zero  set  of  Z{f).  Consider  the  quotients  (11).  The  smoother  the 
quotient,  the  more  rapid  the  decay  of  the  Gabor  coefficients  in  the  sense  described 
by  harmonic  analysis  theorems. 

Computer  experiments  have  verified  these  comments  in  the  special  case  of  the 
Gaussian  g  with  signals  /  taken  from  the  Hermite  functions.  It  should  be  pointed 
out  that  although  the  theory  is  about  analog  signals  computation  must  take  place  on 
the  discrete  level. 

If  the  Zak  transform  of  /  has  two  zeros  in  the  unit  square  then  the  harmonic 
analysis  criteria  can  point,  say,  to  a  Gaussian  type  g  having  exactly  one  of  these  zeroes 
in  the  unit  square  but  from  a  more  pattern  recognition  point  of  view  it  might  better  to 
choose  a  g  whose  Zak  transform  zeroes  match  exactly  the  zero  set  of  Z{f).  Preliminary 
studies  bear  out  this  possibility.  A  Hermite  function  having  exactly  three  zeroes  in 
the  unit  square  was  taken  as  the  the  analyzing  signal.  Higher  order  Hermite  functions 
were  expanded  relative  to  the  Gaussian  and  relative  to  the  analyzing  Hermite  function. 
The  Hermite  analyzing  signal  produced  a  more  compact  representation  in  most  cases. 

When  available,  Tayk  r  series  expansions  at  the  zeroes  can  be  matched  as  much 
as  possible. 

A. 3  Generalized  Gabor  Expansions 

A  signal  /  will  not  always  have  a  Gabor  expansion  relative  to  an  analyzing  signal 
g.  The  selection  procedure  of  the  previous  section  establishes  rules  for  degrees  of 
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compatibility  between  /  and  g  based  on  the  comparison  of  their  Zak  transforms, 
either  by  smoothness  of  the  quotient  or  optimal  matching  of  zeroes.  However,  it  can 
happen  in  some  apphcations  that  a  fixed  g  must  be  taken.  The  problem  is  then  to 
modify  the  definition  of  Gabor  expansions  so  that  an  /  can  be  expanded  as  a  modified 
Gabor  expansion.  Again  our  goal  is  to  produce  an  adaptive  procedure  which  depends 
on  both  /  and  g. 

The  main  idea  is  that  the  definition  given  of  Gabor  expansion  implicitly  depends 
on  the  lattice  Z  and  that  by  extending  the  lattice,  we  can  guarantee  Gabor  expansions. 
More  precise  results  for  the  Gaussian  case  appear  on  [30]  which  include  algorithms  of 
Gabor  coefficient  computation  will  be  described  below. 

Set  5^1  =  e~‘^^  and  g2  =  .  The  Zak  transforms  of  gi  and  <72 

vanish  uniquely  at  the  point  x  =  y  =  ^  and  x  =  j/  =  0,  respectively,  in  the  unit 
square.  Denote  the  r-times  continuously  differentiable  functions  in  the  plane  by  C”". 
The  first  result  we  have  is  that  if  Z{f)  €  C"',  then 

Z{f)  =  pZ{gi)  +  qZ{g2), 

where  p  and  q  are  trigonometric  series  in  C"',  i.e., 

?(!,!/)=  E  (12) 

m^Z  ti€^Z 

i7i€  Z  n€  Z 

If  r  <  2,  say,  then  we  can  write 

/=EE  ^mn  (^1)  mn  EE  ^mn(^92)mn^  (13) 

m^Z n^Z  mgZ  neZ 

where 

I  «mn  r<  00, 

m  n 

I  P<  00. 

m  n 

Since  g2  =  we  can  write  (  13)  as 

f  ~  ^2  ^2  ^rnn9m,^'> 

TFi^Z  n^Z 

where  Cm^in  =  o,mn  and  Cjn2n+  ^mn  The  lattice  of  the  expansion  is  Z  x  jZ.  The 
overall  effect  is  to  double  th-  samp!  ig  rate.  Similar  results  can  be  found  in  the 
language  of  frames  in  [18]. 
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The  proof  of  (  12)  as  given  in  [5]  depends  on  a  partition  of  unity  argument  adapted 
to  the  Heisenberg  group  and  easily  generalizes.  In  particular,  the  choice  of  the  Gaus¬ 
sian  g\  and  g2  are  arbitrary  and  can  be  replaced  by  any  two  distinct  time-frequency 
translates  of  the  Gaussian  g  =  .  In  general,  results  of  the  form  (  12)  depend 

solely  on  the  disjointness  of  the  zero  sets  of  g\  and  g2  and  can  be  extended  to  any 
number  of  analyzing  signals.  Results  of  the  form  (  13)  require  that  these  signals  be 
related  by  time-frequency  translates. 

The  central  result  in  [30]  is  that  q  can  always  be  taken  as  a  trigonometric  poly¬ 
nomial  whose  coefficients  are  explicitly  computable.  Uniqueness  and  algorithms  con¬ 
structing  p  and  q  are  given  in  [30].  A  quantitative  relationship  is  established  between 
the  desired  smoothness  of  p  and  the  degree  of  the  trigonometric  polynomial  q.  Com¬ 
puter  experiments  have  been  carried  out  for  the  Gaussian.  Other  important  special 
functions  whose  role  is  more  than  an  analyzing  signal  but  is  also  part  of  the  applica¬ 
tion  have  been  studied.  In  particular,  Hermite  functions  are  intrinsically  interesting 
in  many  applications  including  image  coding,  computer  vision  and  human  visual  per¬ 
ception  [22]. 

The  main  result  in  [30]  is  that  for  Z{f)  G  we  can  uniquely  write 
Z{f)  ==  pZ{gi)  +  qZ{g2), 

where  p  G  and  ^  is  a  trigonometric  polynomial  of  degree  3  whose  coefficients  are 
explicitly  computable  in  terms  of  the  partial  derivatives  of  Z(/)  at  x  =  ?/  =  |.  Recent 
results  [20]  have  implemented  the  computation  of  p  and  q  soky  using  Zak  transform 
methods. 

The  Fourier  transform  plays  a  major  part  in  [1]  and  [30].  Functions  are  decom¬ 
posed  into  their  eigen  vector  subspaces  relative  to  the  Fourier  transform  which  on  the 
plane  is  given  by  the  linear  operator 


A.4  Digital  Computations 

Gabor  expansions  must  be  finitized  for  digital  computations.  Recent  efforts  [4,  23,  32] 
have  subjected  Gab  or- expansions  to  the  same  periodization  and  sampling  procedures 
which  underlie  the  classical  Fourier  sampling  theory  and  have  contrasted  results  to 
standard  truncation  and  sampling.  The  main  digital  Gabor  expansion  formula  will 
be  described  below  but  the  main  point  to  be  emphasized  is  that  the  analyzing  signal 
translates  undergo  periodization  creating  overlap.  In  the  Fourier  theory,  the  only 
relevant  periodization  occurs  in  the  expansion  coefficients  since  the  basis  functions 
remain  unchanged  under  periodization.  The  appearance  of  periodized  analyzing  sig¬ 
nals  in  the  digital  Gabor  expansion  case  is  a  reflection  of  the  nonlinearity  of  such 
expansions. 

Suppose  a  signal  /  has  a  Gabor  expansion 
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f  —  ^  ^  'y  ]  ^mn9mn 

m€Z  n^Z 

satisfying  Yl,m^z  Yhn^z  I  Omn  P<  oo-  The  Gabor  coefficients  are  uniquely  determined 
by  the  norm-square  finite  energy  condition.  In  general,  the  system 

{gmn  :  m,ne  Z] 

i:;  not  orthogonal  so  straightforward  computation  of  Gabor  coefficients  is  not  always 
possible.  The  design  of  accurate  and  stable  algorithms  for  computing  Gabor  coef¬ 
ficients  (analysis)  and  for  computing  input  signal  from  Gabor  transform  methods 
useful  signal  processing  tools.  This  step  is  intimately  tied  to  the  form  and  meaning 
assigned  to  digital  Gabor  expansions.  We  will  review  briefly  how  periodization  sam¬ 
pling  procedures  force  such  digital  Gabor  expansions.  Such  expansions  as  contrast  to 
more  standard  forms  highlight  the  tradeoff  between  sampling  rate  and  aliasing  errors. 
For  simplicity,  choose  integers  M  >  0  and  iV  >  0  as  sampling  rate  and  periodization 
interval.  Periodize  /  mod  N. 

Mt)  =  E  /('  +  JN), 

)€Z 

and  sample  at  rate 


Q<m<M,  Q<k<N. 


The  resulting  L-tuple  of  values,  L  =  MN,  is  the  digital  signal  corresponding  to  /. 
The  independent  parameters  M  and  N  are  usually  fixed  by  a  priori  information  in  a 
given  application. 

The  samples  (14)  are  related  to  Zak  transform  samples  by 

Z(/)(p^)=E/(i  +  *)e--. 

By  the  second  fundamental  formula,  these  samples  equal 


a=0  r=0 


N-\  M-1 


(15) 


where 


XI  ^r+r>M,>+s'N- 
r'&Zs'eZ 

A  similar  argument  is  the  basis  for  finitizing  the  Fourier  transform.  However,  since 
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Z(a  )(—  —)  =  Z(a)(—  — 


we  can  rewrite  (  15)  as 

^1  n  m  . 

S  ^  "iu)' 

a=0  r=0  •'’  " 

By  formula  (  14),  we  have  the  desired  finite  Gabor  expansion 
N-i  M-l 

Mj7  +  fc)  =  E  E  Ar.{9rM^  +  t). 

a=0  r=0 


0  <  m  <  M,  0  <  71  <  A^. 


(16) 


The  periodized  analyzing  signal  translates 

777 

(5ra)N(^  +  ^),  0<r<M,0<S<N 

form  a  basis  for  signal  expansion  in  L{2ifL),  L  —  MN,  and  were  introduced  in  [4]. 
Algorithms  for  computing  the  Gabor  coefficients  of  the  finite  Gabor  expansion  (  16) 
were  derived  in  the  work  and  increased  resolution  procedures  were  established. 

The  overlapping  of  the  basis  signals  {grs)N  around  the  sampling  interval  intro¬ 
duces  additional  programming  effort  and  at  times  aliasing  errors  which  must  however 
be  accounted  for.  Code  provided  in  this  proposal  will  compare  the  benefits  of  this 
additional  effort  j  standard  truncation-sampling  approaches  and  make  the  optimal 
decision. 


A. 5  New  Algorithms  for  Computing  Gabor  Coefficients 

Several  algorithms  are  by  now  standard  for  computing  Gabor  coefficients.  A  summary 
of  two  appear  in  [3]  along  with  the  basic  constraints  and  tradeoffs. 

The  biorthogonal  approach  introduced  in  [12]  depends  on  solving  the  equation 

Z{h){x,y)Z*{g){x,y)  =  \.  (17) 

The  function  h  is  called  a  biorthogonal  of  g.  If  Z{g)  is  continuous  then  a  solution 
h  €  T^(3?)  need  not  exist.  Neglecting  this  point  for  the  moment,  if  some  h  satisfying 
(  17)  is  found  then  the  Gabor  coefficients  can  be  computed  by  the  formula 

®mn  /?^mn  (^^) 

which  is  an  immediate  result  of  the  first  fundamental  formula. 

The  problem  of  this  approach  for  continuous  Z{g)  is  that  since  h  need  not  be  in 
T^(3?)  the  computation  (  18)  can  be  difficult  to  carry  out.  We  propose  a  generaliza¬ 
tion  which  has  the  effect  of  localizing  Gabor  coefficient  computation  and  producing 
stable  local  computations  by  ‘delaying’  the  instability  to  regions  removed  from  some 
initialization.  We  will  explore  the  approach  for  the  Gaussian  g  = 

Although  (  17)  cannot  be  solved  for  h  €  T*(3J),  we  can  solve 
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Z{h,){x,y)Z*{g){x,y)  =  l^-e^^'^  (19) 

Z(/i2)(x,y)Z*(5)(x,y)  =  l+e2-^  (20) 


for  hx  and  /12  in  Z^(3f).  In  this  case,  ^2  =  the  Fourier  transform  of  hi  but  this 
will  only  be  the  case  when  g  =  g.  From  the  first  fundamental  formula 


1  9mn  ^  — 

h2i  g-mn 


m  =  n  =  0  or  m  =  0,n  =  1 
otherwise, 

m  =  n  =  0  or  m  =  l,n  =  0 
otherwise, 


The  algorithm  proceeds  by  first  precomputing  hi 


and  h2  and  then  computing 


hra  — f^{.^l'}ra  J',  5  €  Z, 

(21) 

CrO=<  f,{hi)ro>,  r  e  Z. 

(22) 

Formula  (  19)  and  (  20)  imply 

bra  ~  (^ra  d"  Qr,3+1) 

(23) 

Crs  —  O'ra  d"  Qr+1,4- 

(24) 

Assuming  that  ooo  is  known,  the  computations  (  21)  and  (  22)  can  be  placed  (  23)  in 
and  (  24)  to  compute  all  0^^. 

The  computations  are  locally  stable  about  coo  but  by  analyzing  the  impulse  re¬ 
sponse  related  to  the  difference  equations  we  see  that  the  resulting  transfer  function 
is  only  marginally  stable.  The  choice  of  initialization  at  coo  is  arbitrary  and  similar 
results  can  be  derived  about  any  initialization.  Stable  global  algorithms  require  a 
precise  understanding  of  the  branching  that  can  occur  at  points  away  from  an  initial¬ 
ization  and  rules  for  providing  a  new  initialization  at  the  ‘boundaries  of  stability’. 

An  increasing  degree  of  smoothness  of  h  can  be  achieved  by  replacing  1  -|- 
in  (  20)  by  trigonometric  polynomials  having  higher  order  zeroes  at  x  =  y  =  1/2. 
However,  the  increasing  degree  of  smoothness  of  h  is  paid  for  by  an  increasing  com¬ 
plexity  of  the  synthesis  algorithm  including  increasing  initializations  and  instability 
as  measured  by  the  impulse  response  poles  on  the  unit  circle  (equal  to  the  degree 
of  the  poles).  The  quantification  of  the  tradoff  between  increased  smoothness  and 
complexity  /  instability  of  synthesis  algorithm  including  measurement  of  bounds  for 
stability  will  be  studied  in  this  proposal. 
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B  The  Ambiguity  Function 

The  results  of  this  section  were  first  proved  in  a  series  of  papers  [26,  27,  28,  29]  in  the 
early  part  of  the  1960s  and  provided  much  of  the  theoretic  framework  during  the  fol¬ 
lowing  decade  for  radar  signal  synthesis.  Throughout,  we  work  with  the  unsymmetric 
form  of  the  cross-ambiguity  function  of  two  signals  f,g€ 

A(f,9){u,v)^  J  (25) 

There  are  two  important  ways  to  view  cross-ambiguity  functions.  First  we  can  con¬ 
struct  the  parametrized  product 

K{t)  =  f{t)gV-v)  (26) 

and  ^'o^m  A{f^g){u,v)  as  the  Fourier  transform  of  h.u{t)  : 

A{f,g){u,v)  =  J  (27) 

In  signal  processing,  A{f,g){u,v)  is  called  the  short-time  Fourier  transform.  The 
function  g{t)  is  viewed  as  a  window  at  the  origin  and  g{t  —  u)  as  a  sliding  window 
which  for  a  particular  o  is  centered  at  v.  The  parametrized  product  hy{t)  presents  a 
windowed  version  of  /  around  v. 

It  is  easy  to  show  that  A{f.,g){u^v)  €  and 

M(/,s)iii=ii/inur-  (28) 

In  fact 

<  A{fi,gi),A{f2,g2)  >=<  fij2  ><  g2,gi  >  ■  (29) 

Also,  we  can  recover  the  windowed  function  by  the  Fourier  inversion  formula 

-v)=  j  A(f,g)(u,v)e^’“'‘du.  (30) 

We  can  go  one  step  further  by  taking  the  Fourier  transform  on  both  sides  with  respect 
to  the  V  variable  : 

firngnyy”'”  =  J  J  A{J,g)(u,v)e^’^^‘'-”Uudv.  (31) 

This  separation  property  under  the  Fourier  transform  h«?>s  been  shown  to  be  necessary 
and  sufficient  for  a  cross-ambiguity  function.  The  left-hand  side  of  (  31)  has  also  been 
identified  as  a  complex  energy  density  function  for  complex  envelope  /  and  g.  That 
is 


mirnyy^^hit 


(32) 
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can  be  interpreted  as  the  differential  complex  energy  in  a  bin  8y  by  8t  located  at  y 
and  t,  where  the  total  complex  energy  is  simply  the  inner  product 

E,=<f,g>=J  f{i)g'{l)dt  (33) 

Equation  (31)  shows  that  this  measure  of  energy  density  may  be  computed  by  evalu¬ 
ating  the  two-dimensional  Fourier  transform  of  the  cross-ambiguity  function  of  /  and 
g  at  the  point  (t,  y).  In  this  light  it  should  not  be  surprising  to  find  A{f,g)  appear  in 
other  ways  that  relate  to  signal  energy  -  see  the  discussion  of  Weyl-Heisenberg  frames 
in  section  4  of  this  paper. 

A  second  way  of  writing  the  cross-ambiguity  function  is  as  an  inner  product 

9){u,  v)=<  f,  >,  (34) 

where  g^vit)  —  dit  —  Since 

^(t)  =  (35) 

we  have  by  the  Plancheral  theorem  that 
=  <f^9^v> 

=  J  f{t)g{t  - 

=  e-2-“M(/,y)(-t;,u). 

If  we  restrict  (u,u)  to  a  lattice  of  points  {mM,nN),  m,n  €  Z,  the  corresponding 
set  of  samples 

A{fi9){'<^M,nN)=<f,gmM,nN>-,  m,n£Z  (36) 

is  called  the  discrete  short-time  Fourier  transform  of  /  relative  to  the  analyzing  signal 
g  and  the  lattice  determined  by  {M,N).  The  set  of  functions 

{5mM,n7V  :  m,n  €  Z},  (37) 

is  usually  called  a  Weyl-Heisenberg  wavelet  system  but  in  this  work  we  will  call  it 
simply  a  Gabor  wavelet  system.  It  will  be  denoted  by  (y,  Af,  N).  We  can  ask  whether 
from  the  set  of  inner  products  (  36)  a  numerically  stable  algorithm  reconstructs  or 
approximates  /.  As  explained  in  [17]  the  question  can  be  answered  affirmatively  if 
and  only  if  the  following  condition  is  satisfied  : 

Frame  Condition  There  exist  constants  0  <  A  <  5  <  oo  such  that  for  all  /  € 

P(X), 

^  II  /  ll'<  E  E  l<  >p<  B  II  /  f  .  (38) 

m  n 

In  this  case,  we  call  {g,M,N)  a  frame  and  A  and  B  frame  bounds  If,  for  some 
constant  0  <  .B  <  oo. 
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EE  l< >I'=B  II /IP  (39) 

m  n 

whenever  /  G  we  say  that  {g,M,N)  is  a  tight  frame.  Condition  (  38)  easily 

translates  into  a  statement  about  cross-ambiguity  function  samples  : 

A  II  /  |P<  EE  I  A(f,g)(mM,nN)  \-‘<  B  ||  /  ||^  .  (40) 

m  n 

Our  tools  for  analyzing  frame  conditions  are  based  on  the  following  analog  sum¬ 
marized  in  [16]. 

Radar  Theorem  (1960) 

Theorem  5  If  f ,  g  €  then 

j /  I  A{!,g)(u,v)  1^  (41) 

=  A{f){y,-x)A*{g){y,-x). 

Theorem  6  If  fi,  /a,  gi,  g2  G  then 

J  j  AUuh){u,v)A-(gug2)Me-^'^^'‘*'''’Uudv  (42) 

=  A{fi,gx){y,-x)A*{f2,g2){y,-x). 


In  [24],  the  Poisson  summation  formula  was  used  to  derive  discrete  analog  of 
theorems  1  and  2.  We  describe  these  discrete  analogs.  The  exact  conditions  on  /  and 
g  required  for  these  formulas  are  given  in  [24]. 


TT^ 

Theorem  7  y^/(t  -F  nN)g*(t  A  nN  —  — ) 

*1  i»J 


1 


=  ^E4l(/,j)(^.^)e^''»‘.  rn€M. 

Theorems  ^  |  44(/,  5^)(mM,nA'’) 

m  n 

=  mEEA(f)(^,^)A-ig)i^,^). 

Theorem  9  YU2Mfu9i){'mM,nN)A*{f2,g2){mM,nN) 

m  n 

=  -^Y.Y,A(fuh)('^,^)A-(g„g,){l,^). 


N'  M‘ 
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C  Applications 

We  say  {g,  jj,  IM)  satisfies  Condition  A,  if  A{g){^,mM)  =  0,  unless  m  =  n  =  0, 
First  we  show  that  there  exists  g  €  satisfying  condition  A.  Define 


9o{t)  =  ^ 

Since  N  < 


1,  0<t<N, 
0,  otherwise. 


m . 


9o{i)g*o{i  “  unless  m  =  0. 

Condition  A  follows  from 

AM(j,0)  =  I" 

N,  n  =  0, 

O,  otherwise. 


('6) 


(47) 


=  ( 


In  [24]  the  following  result  was  proved. 

Theorem  10  {g,M,N)  is  a  tight  frame  if  and  only  if  g  satisfies  condition  A.  In  this 
the  frame  bound  is 


B  = 


MN 


(48) 


The  proo!  of  the  theorem  and  other  results  found  in  [24]  rest  on  the  application  of 
the  discrete  formulas  to  control  the  middle  sum  of  the  frame  condition  by  imposing 
constraints  on  the  values  of  the  auto-ambiguity  of  the  analyzing  signal  g  over  the  dual 
lattice.  Usually,  no  condition  except  for  square- intcgrability  is  imposed  on  /.  The 
application  of  these  ideas  to  the  design  of  adaptive  methods  depends  on  imposing 
joint  conditions  on  the  auto-ambiguities  of  /  and  g  on  the  dual  lattice. 
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